{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LEGACY EXAMPLE: THIS NOTEBOOK IS NO LONGER MAINTAINED AND IS NOT EXECUTABLE WITH UP TO DATE VERSIONS OF THE PACKAGES IT DEPENDS ON\n",
    "\n",
    "# _*Quantum chemistry with Qiskit `Terra`, `Aqua` and CQC's `t|ket〉` compiler*_  \n",
    "\n",
    "In this tutorial, we discuss how to use IBM's Qiskit `Terra` and `Aqua` packages, and the `t|ket〉` compiler by Cambridge Quantum Computing (CQC), to calculate the excited state energies of simple molecules using the quantum subspace expansion (QSE) technique. By the end of this tutorial, you will\n",
    "\n",
    "* Understand why optimizing circuits for quantum chemistry is necessary\n",
    "* Learn about the quantum subspace expansion technique for computing excited state energies\n",
    "* See how Qiskit `Terra` enables 3rd-party passes in its transpiler architecture\n",
    "* Use `pytket` to perform native circuit optimization and circuit routing\n",
    "\n",
    "NOTE: Throughout this tutorial, we assume the reader has some familiarity with quantum chemistry methods, though an extensive knowledge is not necessary.\n",
    "\n",
    "**Code setup**\n",
    "\n",
    "This tutorial makes use of Qiskit `Terra` and `Aqua`, as well as `pytket` (the Python interface to `t|ket〉`). To install `Terra` and `Aqua`, follow instructions available [here](https://qiskit.org/). To install `pytket`, follow the instructions at [this Github repository](https://github.com/CQCL/pytket)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum chemistry in the NISQ era, and the need for circuit optimization techniques\n",
    "\n",
    "One of the expected applications of noisy, intemediate-scale quantum (NISQ) devices is simulating quantum chemistry. A simple quantum chemistry question is \"Given how the electrons and nuclei of a particular molecule are arranged, what is the _energy_ of that configuration?\". By examining how the energy of the configuration changes as the nuclei and electrons are moved relative to one another, we can map out an _energy surface_. The configuration that minimizes the energy is a stable, equilibrium point for the molecule. Knowing this configuration (and its associated energy), we can deduce a variety of molecular properties, such as reaction rates.\n",
    "\n",
    "Most techniques for computing molecular energies on near term quantum devices rely on the [variational quantum eigensolver (VQE) algorithm](https://doi.org/10.1038/ncomms5213). This algorithm uses a parameterized _ansatz_ $|{\\psi}(\\boldsymbol{\\theta})\\rangle$ to describe the ground state energy of the Hamiltonian $H$ for a given configuration. By examining how the expected energy $\\langle \\psi(\\boldsymbol{\\theta})|H| \\psi(\\boldsymbol{\\theta})\\rangle/\\langle \\psi(\\boldsymbol{\\theta})|\\psi(\\boldsymbol{\\theta})\\rangle$ changes as $\\boldsymbol{\\theta}$ is varied, we can optimize the parameters to find an estimate for the ground state and its corresponding energy.\n",
    "\n",
    "Running the VQE algorithm on actual hardware has two complications:\n",
    "\n",
    "* Efficiently preparing the trial state $|\\psi(\\boldsymbol{\\theta})\\rangle$.\n",
    "\n",
    "* Efficiently measuring the terms in the Hamiltonian $H$ that describe the configuration.\n",
    "\n",
    "In this tutorial, we'll focus on the problem of efficiently preparing the trial state. We do so for two reasons:\n",
    "\n",
    "* An accurate estimate of the ground state is all that's necessary to do the quantum subspace expansion technique to estimate excited state energies\n",
    "\n",
    "* To demonstrate the capabilies of Qiskit `Terra` and CQC's `t|ket〉`  compiler to reduce the circuit _resources_ (number of gates, depth, etc.) for the state preparation. In principle, these capabilities can be deployed for other problems.\n",
    "\n",
    "Reducing the resources is crucial in the NISQ era, because the noise present in NISQ devices dominates more for larger circuits (with more operations) and reduces the accuracy of the results. Therefore techniques which reduce circuit requirements can significantly improve the results on real hardware. Note that there are techniques, such as [qubit tapering](https://arxiv.org/abs/1701.08213), or ['gate-efficient' circuits](https://arxiv.org/abs/1809.05057), which are useful for constructing a lower-resource circuit that could be run on _ideal_ hardware. In practice, even these \"resource-efficient\" circuits can be further optimized, especially when the imperfections and constraints of a real piece of hardware are taken into account.\n",
    "\n",
    "At a high level, our approach utilizes Qiskit `Aqua` to generate a parameterized ansatz based on the [_Unitary Coupled Cluster, Single-Double_ (UCCSD) ansatz](https://en.wikipedia.org/wiki/Coupled_cluster). We then use Qiskit `Terra` and the `t|ket〉` compiler (as implemented in `pytket`) to optimize the circuit for that ansatz. We will show that the optimized ansatz requires substantially fewer circuit resources than the naive UCCSD ansatz. An application of our approach, we show how to use the quantum subspace expansion (QSE) technique to compute excited state energies for gaseous hydrogen (H$_{2}$) and lithium hydride (LiH)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing excited state energies using the Quantum Subspace Expansion (QSE) technique\n",
    "\n",
    "The VQE algorithm is often used to compute the _ground state_ energy of a given molecular configuration. However, knowing the energies of _excited states_ of the system is also useful. The energies of electronic excited states prove, in general, challenging to compute. Excited states are usually more _entangled_ than ground states and so require more computational resources to compute. While there exists techniques for efficiently representing certain classes of entangled states, in general, the complexity of the simulation will be non-trivial. This is a problem as it makes classical computation of the energy, to an appropriate accuracy, of many molecules impossible. \n",
    "\n",
    "This notebook demonstrates calculation of excited states of molecules using the [quantum subspace expanansion (QSE) technique](https://doi.org/10.1103/PhysRevA.95.042308). The QSE technique uses an accurate estimate of the ground state energy of a given molecular configuration to estimate the energies of excited states. Consider a fixed molecule (represented by a given Hamiltonian $H$), and suppose $\\left|\\Psi_{0}\\right\\rangle$ is the output of the VQE algorithm for estimating the ground state energy of $H$.\n",
    "\n",
    "The QSE technique constructs a subspace of state vectors $\\left|\\Psi_j^k\\right\\rangle$ formed by one-electron excitations of the ground state wavefunction:\n",
    "\n",
    "\\begin{equation}\n",
    "\\left|\\Psi_{j}^{k}\\right\\rangle = c_k^{\\dagger}c_{j}\\left|\\Psi_0\\right\\rangle.\n",
    "\\end{equation}\n",
    "\n",
    "where $c_k^{\\dagger}, c_{j}$ are the fermionic creation and annihilation operators over spin orbitals $k$ and $j$, respectively. That is, these vectors are formed by reducing the occupation of spin orbital $j$ by one, and increasing the occupation of spin orbital $k$ by one. The vectors are not in general orthogonal to $\\Psi_{0}$ hence we will need to calculate an overlap matrix.\n",
    "\n",
    "\n",
    "Within this subspace, we solve a generalized eigenvalue problem. Consider the operator $H'$ with matrix elements given by\n",
    "$$(H')_{jk}^{lm} = \\langle\\Psi_j^l \\left| H \\right| \\Psi_k^m\\rangle,$$\n",
    "\n",
    "and define an overlap matrix $S$ whose matrix elements are given by\n",
    "\n",
    "$$S_{jk}^{lm} = \\langle \\Psi_j^l \\left|\\Psi_k^m\\right\\rangle.$$\n",
    "\n",
    "The generalized eigenvalue equation to be solved is \n",
    "\n",
    "\\begin{equation}\n",
    "H'C=SCE,\n",
    "\\end{equation}\n",
    "\n",
    "where $C$ is the matrix of eigenvectors, and $E$ is the vector of eigenvalues. Crucially, _the energy eigenvalues $E$ provide an estimate of the excited state energies of $H$_ as well as a refined value of the ground state energy.\n",
    "\n",
    "Notice that the solution to the generalized eigenvalue equation can be done on a classical computer, provided $H'$ and $S$ have been calculated. The matrix elements of both of these matrices can be constructed using a quantum computer, in the following way. First, re-write the matrix elements in terms of $\\left | \\Psi_{0}\\right \\rangle$:\n",
    "\n",
    "\\begin{align}\n",
    "(H')_{jk}^{lm} =& \\langle\\Psi_j^l \\left| H \\right| \\Psi_k^m\\rangle = \\langle \\Psi_{0} | c_{j}^\\dagger c_{l} Hc_{m}^{\\dagger}c_{k}|\\Psi_{0}\\rangle\\\\\n",
    "S_{jk}^{lm} &= \\langle \\Psi_j^l \\left|\\Psi_k^m\\right\\rangle = \\langle \\Psi_{0} | c_{j}^\\dagger c_{l} c_{m}^{\\dagger}c_{k}|\\Psi_{0}\\rangle.\n",
    "\\end{align}\n",
    "\n",
    "The matrix elements can be calculated using a quantum computer or simluator. How? By transforming the operators \n",
    "$c_{j}^\\dagger c_{l} c_{k}^{\\dagger}c_{m}$ and $c_{l}^\\dagger c_{j} H c_{m}^{\\dagger}c_{k}$ to a set of Pauli quantum gates according to an appropriate scheme such as Jordan-Wigner or Bravyi-Kitaev, apply this gate set to the ground state wavefunction (constructed with the coefficients obtained from the VQE calculation) and perform a measurement to obtain the expected value. By measuring these expectation values, we obtain estimates of $S_{jk}^{lm}$ and $(H')_{jk}^{lm}$, respectively.\n",
    "\n",
    "Therefore, to use the QSE technique, we need to use the quantum computer to calculate 2 quantities:\n",
    "\n",
    "* An estimate of the ground state $\\left |\\Psi_{0}\\right\\rangle$ (by running the VQE algorithm)\n",
    "* The matrix elements of $H'$ and $S$.\n",
    "\n",
    "For the first quantity, we want to have an efficient representation of the trial state $|\\psi(\\boldsymbol{\\theta})\\rangle$. And for the second, we need to be able to take the estimate of the ground state and efficiently compute matrix elements. In both of these cases, being able to optimize the circuit for preparing would be useful. Thankfully, Qiskit `Terra`, in conjunction with `t|ket〉`, provides us tools for doing so. We discuss the problem of _circuit compilation/optimization_ in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circuit compilation with Qiskit `Terra` and `t|ket〉`\n",
    "\n",
    "In the remainder of this tutorial, we show how to use Qiskit `Terra` and `t|ket〉` to optimize the circuits for preparing VQE trial states, and simulate the QSE technique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application I: Reducing circuit resources for trial state preparation\n",
    "\n",
    "In this tutorial, we'll focus on two simple molecules: hydrogen (H$_{2}$) and lithium hydride (LiH). NOTE: the code is much slower for LiH, which is why here, we'll demonstrate H$_{2}$.\n",
    "\n",
    "\n",
    "As a first application of the pipelines provided by `Terra`, `Aqua`, and `t|ket〉`, we first show how to configure a VQE experiment using `Aqua`, and then use `Terra` and `t|ket〉` to compile the trial state preparation circuit.\n",
    "\n",
    "Let's start by importing the necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# numpy, for random number generation\n",
    "import numpy as np\n",
    "\n",
    "# Qiskit, for transpiler-related functions, the IBMQ provider, and the Aer simulator\n",
    "from qiskit import IBMQ, Aer, QuantumRegister\n",
    "from qiskit.transpiler import transpile, transpile_dag, PassManager\n",
    "from qiskit.converters import circuit_to_dag, dag_to_circuit\n",
    "\n",
    "# pytket, for optimization\n",
    "import pytket\n",
    "from pytket.qiskit import TketPass\n",
    "\n",
    "# Qiskit Aqua, for chemistry\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType\n",
    "from qiskit.chemistry import FermionicOperator\n",
    "from qiskit.chemistry.components.initial_states import HartreeFock\n",
    "from qiskit.chemistry.components.variational_forms import UCCSD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 0: Enable IBMQ account"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-29T00:04:16.313210Z",
     "start_time": "2018-09-29T00:04:14.460647Z"
    }
   },
   "outputs": [],
   "source": [
    "provider = IBMQ.load_account()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: setting up the molecule\n",
    "\n",
    "We choose the basis set spanning the molecular wavefunction, the molecular geometry, the chemical identity of each atom, the charge and spin quantum number.\n",
    "\n",
    "To calculate results for LiH (slower), just comment out the H$_{2}$ string and replace it with the LiH one.\n",
    "\n",
    "NOTE: Here, we focus only on one particular value for the `bond_length`. If you wanted to replicate the final plot for the excited state energies of LiH as a function of bond length, you'd need to sweep `bond_length` over several values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-29T00:04:20.069592Z",
     "start_time": "2018-09-29T00:04:20.065489Z"
    }
   },
   "outputs": [],
   "source": [
    "# Choose a particular bond length\n",
    "# NOTE: Units are in Angstroms\n",
    "\n",
    "bond_length = 0.7\n",
    "\n",
    "# Set up molecule\n",
    "\n",
    "# base_molecule_str = 'Li .0 .0 .0; H .0 .0 {}'\n",
    "base_molecule_str = 'H .0 .0 .0; H .0 .0 {}'\n",
    "\n",
    "# Specify other molecular properties\n",
    "charge = 0\n",
    "spin = 0\n",
    "basis = 'sto3g'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having set up the molecule, we now execute our classical chemistry driver to obtain the integrals that define the terms in the molecule's Hamiltonian. In this case, we choose `PYSCF` as our driver, so make sure you have that installed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Molecular repulsion energy:  0.7559674441714287\n"
     ]
    }
   ],
   "source": [
    "# Using driver to get fermionic Hamiltonian\n",
    "# PySCF example\n",
    "\n",
    "driver =  PySCFDriver(atom=base_molecule_str.format(bond_length),\n",
    "                      unit=UnitsType.ANGSTROM,\n",
    "                      charge=charge,\n",
    "                      spin=spin,\n",
    "                      basis=basis)\n",
    "\n",
    "molecule = driver.run()\n",
    "\n",
    "print(\"Molecular repulsion energy: \", molecule.nuclear_repulsion_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The molecular repulsion energy calculated by the driver corresponds to the coulombic repulsion between the nuclei in the molecule, and we can add it to our electron structure calculation at the end to get the total energy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Set up a variational form for VQE\n",
    "\n",
    "To run the VQE algorithm, we need to specify a mapping from the molecular Hamiltonian to qubits, an initial state, and the ansatz we use for the trial state. Here, we use the Jordan-Wigner transform to map the molecular Hamiltonian onto qubits (Pauli operators). We choose the initial state to be a Hartree-Fock state, and we take the variational ansatz to be the unitary coupled cluster with single and double excitations (UCCSD)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of qubits: 4\n"
     ]
    }
   ],
   "source": [
    "n_qubits = molecule.one_body_integrals.shape[0]\n",
    "n_electrons = molecule.num_alpha + molecule.num_beta - molecule.molecular_charge\n",
    "\n",
    "# get fermionic operator and mapping to qubit operator\n",
    "ferOp = FermionicOperator(h1=molecule.one_body_integrals, h2=molecule.two_body_integrals)\n",
    "\n",
    "qubitOp = ferOp.mapping(map_type='JORDAN_WIGNER', threshold=0.00000001)\n",
    "qubitOp.chop(10**-10)\n",
    "\n",
    "# Instantiate the initial state as a Hartree-Fock state\n",
    "initial_hf = HartreeFock(num_qubits=n_qubits, num_orbitals=n_qubits, \n",
    "                    qubit_mapping='jordan_wigner', two_qubit_reduction=False, num_particles= n_electrons)\n",
    "\n",
    "# Create the variational form\n",
    "var_form = UCCSD(num_qubits=n_qubits, num_orbitals=n_qubits, \n",
    "                num_particles=n_electrons, depth=1, initial_state=initial_hf, qubit_mapping='jordan_wigner')\n",
    "\n",
    "# How many qubits do we need?\n",
    "print('Number of qubits: {0}'.format(n_qubits))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Transpile circuit and examine circuit properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, we have instantiated a variational form for the molecule according to the UCCSD formulation. Aqua's `UCCSD` method has returned back to us an abstraction of the VQE variational form. Let's query some of its properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, (-3.141592653589793, 3.141592653589793))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Query the variational form for the number of parameters, and the parameter bounds.\n",
    "var_form.num_parameters, var_form.parameter_bounds[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a VQE circuit, `var_form` has some number of parameters, and each parameter can take values in $[-\\pi, \\pi]$.\n",
    "These parameters control the angles of rotation in the quantum circuit that represents the variational anatz.\n",
    "\n",
    "For a particular set of parameters, there is an assocated quantum circuit. Let's input some fiducial parameter values and query properties of the resulting circuit to introduce some nomenclature for describing circuits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'size': 150,\n",
       " 'depth': 83,\n",
       " 'width': 4,\n",
       " 'bits': 0,\n",
       " 'factors': 1,\n",
       " 'operations': {'u3': 42, 'u2': 40, 'cx': 56, 'u1': 12}}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Instantiate a concrete instance of the VQE ansatz by setting all the parameters to the\n",
    "# arbitrarily-chosen value of 0.\n",
    "var_circ = var_form.construct_circuit(np.zeros(var_form.num_parameters))\n",
    "\n",
    "# Use Terra to convert the circuit to its directed, acyclic graph (DAG) representation.\n",
    "var_circ_dag = circuit_to_dag(var_circ)\n",
    "\n",
    "# The .properties() method of the DAG to get circuit properties.\n",
    "var_circ_dag.properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These circuit properties are:\n",
    "\n",
    "\n",
    "* `size`: The total number of gates in the circuit\n",
    "* `depth`: The total number of _layers_ in the circuit\n",
    "* `width`: The number of qubits in the circuit\n",
    "* `bits`: The number of classical bits in the circuit. (NOTE: Because the circuit prepares a VQE trial state, and does not have any measurements, `bits` will be 0.)\n",
    "* `factors`: The number of tensor factors the circuit could be decomposed into (by looking at the number of weakly connected components of the DAG)\n",
    "\n",
    "The `.properties()` method of the DAG representation also breaks down the total number of gates (`size`) by the individual gates themselves. These are the $u1,u2,u3$ and CNOT gates described [here](https://qiskit.org/documentation/terra/summary_of_quantum_operations.html). You can verify that the total number of gates is in fact equal to `size`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Transpile the circuit using Terra and `t|ket〉`\n",
    "\n",
    "Having instantiated a VQE variational form and examined some of its properties, we'd like to optimize those properties so that the circuit could be run on a near-term device. Terra provides a framework (the [_transpiler_](https://qiskit.org/documentation/terra/overview.html#transpiler)) for manipulating circuits according to certain _passes_, and where the execution of the passes is orchestrated by a _PassManager_, which we can use to do this optimization. Importantly, _the transpiler manipulates the circuit to change circuit properties, without actually changing the input-output relationship the circuit defines_. That is, the transpiler takes a circuit and re-writes it, but doesn't change what the circuit actually does.\n",
    "\n",
    "CQC has written several passes for manipulating quantum circuits which are available via `pytket`. For an extensive discussion of the framework `t|ket〉` uses to manipulate circuits, see **[TODO: include link]**. Here, we'll demonstrate using Terra and `pytket` to optimize a randomly-chosen realization of the VQE variational form. Currently, passes in the transpiler require a backend in order to run. For simplicity, we start by using a simulator backend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab an Aer backend\n",
    "aer_backend = Aer.get_backend('qasm_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose a random set of parameters\n",
    "seed = 0\n",
    "np.random.seed(seed)\n",
    "params = np.random.uniform(low=-3.1, high=3.1, size=var_form.num_parameters)\n",
    "\n",
    "# Construct a random instance of the variational circuit\n",
    "var_circuit = var_form.construct_circuit(params)\n",
    "\n",
    "# Turn the circuit into a DAG\n",
    "var_dag = circuit_to_dag(var_circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This randomly-chosen realization of the variational form has the same circuit properties as the circuit we instantiated in Step 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'size': 150,\n",
       " 'depth': 83,\n",
       " 'width': 4,\n",
       " 'bits': 0,\n",
       " 'factors': 1,\n",
       " 'operations': {'u3': 42, 'u2': 40, 'cx': 56, 'u1': 12}}"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var_dag.properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we set up a transpiler using Terra and the `TketPass` from `pytket`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Terra PassManager object\n",
    "tk_pass_manager = PassManager()\n",
    "\n",
    "# Set up the TketPass\n",
    "tk_pass = TketPass(aer_backend)\n",
    "\n",
    "# Add the TketPass to the PassManager\n",
    "tk_pass_manager.append(tk_pass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the transpiler set up and the realization of the variational form put into a DAG, we can now use the `transpile_dag` function to run the PassManger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_dag_transpiled = transpile_dag(var_dag, pass_manager=tk_pass_manager)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the properties of this circuit to see how the transpiled circuit differs from the original one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'size': 95,\n",
       " 'depth': 57,\n",
       " 'width': 4,\n",
       " 'bits': 0,\n",
       " 'factors': 1,\n",
       " 'operations': {'u3': 30, 'cx': 52, 'u1': 12, 'u2': 1}}"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var_dag_transpiled.properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `t|ket〉` compiler has optimized the circuit, and as we might hope for, _the transpiled circuit has a lower `size` and `depth` than the original circuit._\n",
    "\n",
    "\n",
    "The table below gives properties of the transpiled circuit for H$_{2}$ and LiH (with `seed=0`, corresponding to the worst-case and also most common performance of the transpiler).\n",
    "\n",
    "\n",
    "|      Molecule:  H$_2$         | Total Gates | Depth |  CNOT Count |\n",
    "|--------------------------------------|-------------|---------------|--------------------|\n",
    "| Input circuit                        | 150       | 83          | 56              |\n",
    "| `tket` circuit optimization   (Aer backend)       | 95        | 57          | 52              |\n",
    "\n",
    "\n",
    "|      Molecule:  LiH           | Total Gates |  Depth |  CNOT Count |\n",
    "|--------------------------------------|-------------|---------------|--------------------|\n",
    "| Input circuit                        | 13700       | 9342          | 8064               |\n",
    "| `tket` circuit optimization (Aer backend)         | 7411        | 4416          | 5096               |\n",
    "\n",
    "For both H$_2$ and LiH, circuit optimization using `t|ket〉` reduces the number of gates necessary to prepare the trial state. The depth also decreases, which makes running the circuit more feasible on near-term devices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: Route the circuit onto real hardware"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though we've optimized the circuit, there's no guarantee that it can run, as written, on a real backend. This is because the directionaly of the CNOT gates on the backend may not be respected by the circuit. For this reason, we need to re-write the circuit in such as way that the CNOT gates respect the _coupling map_ of the backend. `t|ket〉` knows how to do this, and handles this problem (called circuit \"routing\") when a real backend is put into the `TketPass` object. Routing refers to the process of making quantum circuits hardware compliant by the addition of SWAP gates such that all multi-qubit interactions occur on adjacent physical qubits.\n",
    "\n",
    "Because we know _a priori_ that the H$_{2}$ and LiH molecules requires at most 12 qubits, we make sure to use an IBMQ backend with no less than 12 qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<IBMQBackend('ibmq_16_melbourne') from IBMQ()>,\n",
       " <IBMQBackend('ibmq_qasm_simulator') from IBMQ()>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Grab only backends that have at least 12 qubits\n",
    "provider.backends(filters=lambda x: x.configuration().n_qubits >= 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the `ibmq_16_melbourne` backend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "real_backend = provider.get_backend('ibmq_16_melbourne')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To route the variational circuit onto real hardware, we simply need to change the backend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Terra PassManager object\n",
    "tk_pass_manager = PassManager()\n",
    "\n",
    "# Set up the TketPass\n",
    "tk_pass = TketPass(real_backend)\n",
    "\n",
    "# Add the TketPass to the PassManager\n",
    "tk_pass_manager.append(tk_pass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can transpile the DAG representation of the variational circuit to a backend-compliant circuit using the `tranpsile_dag` function. Here, because the backend has a non-trivial coupling map, the `TketPass` will perform both circuit optimization and optimal routing calculations. We need to add a register containing the ancilla qubits on the architecture that `TketPass` can use in routing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "blank_qubits = QuantumRegister(len(real_backend.properties().qubits) - var_dag.width())\n",
    "var_dag.add_qreg(blank_qubits)\n",
    "var_dag_transpiled = transpile_dag(var_dag, pass_manager=tk_pass_manager)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'size': 108,\n",
       " 'depth': 69,\n",
       " 'width': 4,\n",
       " 'bits': 0,\n",
       " 'factors': 1,\n",
       " 'operations': {'u3': 30, 'cx': 52, 'u1': 12, 'u2': 14}}"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var_dag_transpiled.properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare these results to Qiskit's own default transpilation by passing in the backend's `coupling_map` to the `transpile_dag` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'size': 119,\n",
       " 'depth': 89,\n",
       " 'width': 14,\n",
       " 'bits': 0,\n",
       " 'factors': 11,\n",
       " 'operations': {'u3': 15, 'cx': 56, 'u1': 14, 'u2': 34}}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transpile_dag(var_dag, coupling_map=real_backend.configuration().coupling_map).properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table below shows circuit properties for H$_{2}$ (with `seed=0`).\n",
    "\n",
    "|      Molecule:  H$_2$         | Total Gates | Overall Depth | Overall CNOT Count |\n",
    "|--------------------------------------|-------------|---------------|--------------------|\n",
    "| Input circuit                        | 150       | 83          | 56              |\n",
    "| `tket` circuit optimization (Aer backend)         | 95        | 57          | 52               |\n",
    "| `Qiskit ` default routing         (real backend)           | 119      | 89         | 56              |\n",
    "| `tket` circuit optimzation + routing (real backend) | 108       | 69        | 52               |\n",
    "\n",
    "The table below shows circuit properties for LiH (with `seed=0`).\n",
    "\n",
    "|      Molecule:  LiH           | Total Gates | Overall Depth | Overall CNOT Count |\n",
    "|--------------------------------------|-------------|---------------|--------------------|\n",
    "| Input circuit                        | 13700       | 9342          | 8064               |\n",
    "| `tket` circuit optimization      (Aer backend)    | 7411        | 4416          | 5096              |\n",
    "| `Qiskit ` default routing           (real backend)         | 42178       | 23367         | 17977              |\n",
    "| `tket` circuit optimzation + routing (real backend) | 19256       | 10022         | 8711               |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application II: using QSE to compute excited state energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous application showed us how to use `t|ket〉` and Terra to optimally route circuits onto real hardware. In the introduction, we observed that in order to use the quantum subspace expansion to compute excited state energies, we need to first come up with an estimate of the ground state. For this, we use VQE.\n",
    "\n",
    "In this application, we'll instantiate a VQE circuit, run it, and use the estimated ground state as input to `pytket`'s quantum subspace expansion function(s). NOTE: We'll use a simulator backend, and will not run the VQE algorithm on real hardware."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code imports\n",
    "\n",
    "# From Aqua, we need \n",
    "from qiskit.aqua import QuantumInstance\n",
    "\n",
    "from qiskit.aqua.algorithms.adaptive import VQE\n",
    "from qiskit.aqua.components.optimizers import L_BFGS_B\n",
    "\n",
    "# From pytket, we need QSE functions\n",
    "\n",
    "from pytket.chemistry import QseMatrices, QSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = Aer.get_backend('statevector_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "pass_manager = PassManager()\n",
    "tk_pass = TketPass(backend)\n",
    "pass_manager.append(tk_pass)\n",
    "\n",
    "quantum_instance = QuantumInstance(backend, pass_manager=pass_manager)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Use VQE to estimate ground state\n",
    "\n",
    "First, we'll use VQE to estimate the ground state and its energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Temporary code for Aer on Macbook\n",
    "import os\n",
    "os.environ['KMP_DUPLICATE_LIB_OK']='True'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GS Minimum value: -1.1361894540653963\n",
      "GS Parameters: [ 5.06008657e-07  5.12457730e-07 -1.04867316e-01]\n"
     ]
    }
   ],
   "source": [
    "# Set initial values of parameters\n",
    "number_amplitudes = len(var_form._single_excitations)+ len(var_form._double_excitations)\n",
    "\n",
    "amplitudes_0 = []\n",
    "for i in range(number_amplitudes):\n",
    "    amplitudes_0.append(0.00001)\n",
    "\n",
    "optimizer = L_BFGS_B()\n",
    "optimizer.set_options(maxfun=1000, factr=10, iprint=10)\n",
    "\n",
    "# setup VQE with operator, variation form, and optimzer\n",
    "vqe_algorithm = VQE(operator=qubitOp, \n",
    "                    var_form=var_form, optimizer=optimizer, initial_point=amplitudes_0)\n",
    "\n",
    "results = vqe_algorithm.run(quantum_instance)\n",
    "\n",
    "eigval = results['eigvals'][0]\n",
    "gs_energy = eigval.real + molecule.nuclear_repulsion_energy\n",
    "\n",
    "print(\"GS Minimum value: {}\".format(gs_energy))\n",
    "print(\"GS Parameters: {}\".format(results['opt_params']))\n",
    "\n",
    "# store ground state amplitudes for subsequent steps\n",
    "opti_amplitudes = results['opt_params']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Use QSE to find excited states from ground state <a id='sectionB'></a>\n",
    "\n",
    "We now have the main ingredients to perform a QSE calculation: the molecular Hamiltonian and the optimized parameters to reconstruct the ground state wavefunction. We build our excitation hamiltonian and overlap operators, and measure the elements that compose $H$ and $S$. After we have obtained these arrays, we perform a diagonalization to obtain the excited state energies and vectors.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Excited State Energies:  [-1.13618945 -0.47845306 -0.47845306 -0.47845306 -0.1204519   0.5833141\n",
      "  0.75596744  0.75596744  0.75596744  0.75596744  0.75596744  0.75596744\n",
      "  0.75596744  0.75596744  0.75596744  0.75596744]\n"
     ]
    }
   ],
   "source": [
    "qubitOp = ferOp.mapping(map_type='JORDAN_WIGNER', threshold=0.00000001)\n",
    "n_qubits = qubitOp.num_qubits\n",
    "qubitOp.chop(10**-10)\n",
    "\n",
    "# Use matrix term helper class\n",
    "matrix_terms = QseMatrices(qubitOp, n_qubits)\n",
    "\n",
    "# Instantiate an instance of the QSE algorithm\n",
    "qse_algorithm = QSE(matrix_terms, 'matrix', var_form, opt_init_point=opti_amplitudes)\n",
    "\n",
    "# Run the algorithm\n",
    "energies = qse_algorithm.run(quantum_instance)['eigvals']\n",
    "\n",
    "# The excited state energies are the energies from above,\n",
    "# plus the nuclear repulsion energy.\n",
    "print(\"Excited State Energies: \", energies+molecule.nuclear_repulsion_energy)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculation provides a refined value of the ground state and a series of excited state energies, whose number depends on the size of the basis set chosen for the molecule, as well as its nature and symmetry. Some of the energies obtained are repeated several times, signaling that some of the states obtained are degenerate. This result can be improved by either increasing the basis set or considering higher-order excitations in the subspace expansion.\n",
    "\n",
    "The following graph shows us the excited states of LiH at a range of bond distances calculated via our method, compared to values computed using the classical [EOM-CCSD method](https://aip.scitation.org/doi/10.1063/1.464746). To generate this data yourself, you can  scan the bond length parameter we set at the start of the calculation.\n",
    "\n",
    "![alt text](LiH.png \"Title\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the QSE technique, and our simulated VQE ground state, we find the ground state curve has a minimum at a separation of about 1.5 Å, which is in reasonable agreement with experimental data. The calculation also finds a number of excited states. Looking at the first three of these, we find that at the equilibrium distance, these states are 0.11, 0.12 and 0.17 Ha higher in energy than the ground state, which is again in reasonable agreement with experimental data. Note the small kink in one of the excited states energy curve at a distance of approximately 1.2 Å. This indicates that our restriction to single electron excitations is not enough to provide an accurate description at this distance. Overall, the comparison with classically computed EOM-CCSD curves shows that this method reproduces excited state energies with good accuracy at most distances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}